library(sasld)

# -------------------------------------------------------------------
# ATTENZIONE: la funzione sample e' stata modificata 
# a partire da R 3.6.0 per cui i "vecchi risultati"
# non posso essere riprodotti di default
# http://hostingwin.unitn.it/micciolo/PeM/01-rng.pdf
# -------------------------------------------------------------------
# se si vogliono riprodurre i risultati dei Laboratori
# con le versioni pił recenti di R va eseguita la seguente istruzione
# -------------------------------------------------------------------
RNGkind(sample.kind="Rounding")
# -------------------------------------------------------------------


# 5.1 ---------------------------------------------------------------
dado <- c(1:6)

sample(dado,100)

set.seed(246135)
out <- sample(dado,100,replace=TRUE)

out[10]

sum(out == 6)

out <- sample(dado,100,replace=TRUE)
sum(out == 6)

set.seed(123456)
for (i in 1:20) {
    out <- sample(dado,100,replace=TRUE)
    cat(sum(out == 6)," ")
}

x <- dbinom(c(0:100),100,1/6)
sum(x[24:101])
sum(x[1:13])
sum(x[1:13])/sum(x[24:100])


set.seed(123456)
for (i in 1:20) {
    out <- sample(dado,1000,replace=TRUE)
    cat(sum(out == 6)," ")
}


set.seed(123456)
for (i in 1:20) {
    out <- sample(dado,10000,replace=TRUE)
    cat(sum(out == 6)," ")
}
# -------------------------------------------------------------------




# 5.2 ---------------------------------------------------------------
num <- c(364:341)
den <- rep(365,length(num))
p <- num/den
1 - prod(p)

giorni <- c(1:365)
set.seed(123456)
out <- sample(giorni,25,replace=TRUE)
out

duplicated(out)

x <- c(1:25); y <- out; tmp <- cbind(x,y)
o <- order(y); tmp[o,]


set.seed(654321)
out <- sample(giorni,25*100000,replace=TRUE)
X <- matrix(out,ncol=25)
Y <- apply(X,1,duplicated)
ris <- colSums(Y)
sum(ris > 0)

set.seed(654321)
out <- sample(giorni,23*100000,replace=TRUE)
X <- matrix(out,ncol=23)
Y <- apply(X,1,duplicated)
ris <- colSums(Y)
sum(ris > 0)

set.seed(654321)
out <- sample(giorni,22*100000,replace=TRUE)
X <- matrix(out,ncol=22)
Y <- apply(X,1,duplicated)
ris <- colSums(Y)
sum(ris > 0)
# -------------------------------------------------------------------




# 5.3 ---------------------------------------------------------------
scatola <- c(1,1,1,1,1,1,1,5,10,20)

set.seed(123456)
sample(scatola,2,replace=FALSE)

set.seed(123456)
sum(sample(scatola,2,replace=FALSE))

set.seed(123456)
for (i in 1:20) cat(sum(sample(scatola,2,replace=FALSE))," ")

out <- numeric(100000)
set.seed(654321)
for (i in 1:100000) out[i] <- sum(sample(scatola,2,replace=FALSE))
sum(out >= 20)
# -------------------------------------------------------------------




# 5.4 ---------------------------------------------------------------
gruppo <- c("A","B","0","AB")

set.seed(123456)
sample(gruppo,20,replace=TRUE,prob=c(0.4,0.1,0.45,0.05))

set.seed(123456)
out <- sample(gruppo,20*100000,replace=TRUE,prob=c(0.4,0.1,0.45,0.05))
out <- matrix(out,ncol=20)
job <- function(x) sum(x == "AB")
ris <- apply(out,1,job)
sum(ris == 0)
0.95^20
# -------------------------------------------------------------------